#### Why Do States Admit Refugees? A comparative analysis of resettlement policies in OECD countries

#### LUTZ AND PORTMANN 2021, Journal of Ethnic and migration studies

## Code for replication

## load packages
library(dplyr)
library(xtable)
library(ggplot2)
library(ggpubr)
library(tidyr)
library(grid)
library(gridExtra)
library(sjPlot)
library(estimatr)
library(plm)
library(margins)
library(lmtest)
library(stargazer)
library(plyr)
library(foreign)

## load dataset
load("data_resettlement_oecd.RData")

### Variables:
# pol_openness: resettlement policy openness
# pol_comb: asylum policy mix (liberal resettlement, restrictive control)
# adhoc: adhoc resetttlement program (yes/no)
# quota: resettlement quota (yes/no)
# reset: absolute number of resettled persons
# reset_by_pop: number of resettled persons by 100,000 inhabitants
# pop: population size (in millions)
# gdp_lag: GDP per capita (lagged and logged)
# unemp_lag: unemployment rate (lagged)
# rrpp_vote: radical-right vote share
# lrs: government ideology (left, center, right)
# Refugees: absolute number of refugees worldwide
# refugees_lag: number of refugees worldwide in millions (lagged)
# asylum_seekers: absolute number of asylum seekers
# asylum_requests_lag: number of asylum seekers by 1000 inhabitants (lagged)
# europe --> dummy whether European country or not

# create sub-sample with only European countries
sample_europe <- subset(data_resettlement_oecd, data_resettlement_oecd$europe==1)

# Table A2: Descriptive overview numeric variables
a <- c("pol_openness", "reset_by_pop", "pol_comb", "pop", "gdp_lag", "unemp_lag", "rrpp_vote", "refugees_lag", "asylum_requests_lag")
AT_1 <- data_resettlement_oecd %>% 
  dplyr::select(a) %>% 
  dplyr::mutate_each(funs(as.numeric)) %>%
  summarise_all(list(min = min, mean = mean, max = max, sd = sd), na.rm = TRUE)
AT_1 <- data.frame(matrix(as.vector(AT_1), nrow = 9, ncol = 4))
AT_1 <- sapply(AT_1, as.numeric)
colnames(AT_1) <- c("Min","Mean", "Max", "SD")
rownames(AT_1) <- c("Resettlement policy openness", "Resettlement admissions (share of population size)", "Asylum policy mix", "Population size", "GDP per capita", "Unemployment rate", 
                    "Radical right vote share",  "Refugees worldwide (in millions)", "Asylum request (share of population size)")
AT_1 <- apply(AT_1, c(1,2), function(x) round(x, 2))
xtable(AT_1)

## Figure 1: Comparison of resettlement policies and number of resettled refugees across countries
country_means <- data_resettlement_oecd %>% 
  dplyr::group_by(country) %>% 
  dplyr::summarize(country_m = mean(pol_openness, na.rm=T))
country_means2 <- country_means[order(country_means[,2],decreasing=F),]
resettled_refugees <- data_resettlement_oecd %>% dplyr::group_by(country) %>% 
  dplyr::summarize(reset_by_pop = mean(reset_by_pop, na.rm = TRUE))
resettled_refugees2 <- resettled_refugees[order(resettled_refugees[,2],decreasing=F),]
jpeg("figure1.jpeg", 
     width = 11, height = 5.5, units = "in", res=800)
par (mar=c(3.2,7,1.2,1.2), mgp=c(2,.7,0), tck=-.01, las=1)
par(mfrow=c(1,2))
barplot(country_means2$country_m, horiz=T, main="A. Resettlement policy", xlab="Resettlement policy openness", border=F, names.arg=country_means2$country, las=1, col="black", cex.names=0.7)
abline(v=c(0,0.2,0.4,0.6,0.8,1), lty=2, col="lightgray")
barplot(resettled_refugees2$reset_by_pop, horiz=T, main="B. Resettlement admissions", xlab="Resettlement admissions (per 100,000 inhabitants)", border=F, names.arg=resettled_refugees2$country, las=1, col="black", cex.names=0.7, xlim=c(0,69))
abline(v=c(0,10,20,30,40,50,60), lty=2, col="lightgray")
dev.off()

# Figure 2:
data_resettlement_oecd$pol_openness_f <- factor(data_resettlement_oecd$pol_openness, 
                                                levels=c(0, 0.333333333333333, 0.666666666666667, 1), 
                                                labels=c("no resettlement", "ad-hoc", "quota", "ad-hoc & quota"))
data_resettlement_3periods <- subset(data_resettlement_oecd, data_resettlement_oecd$year==1980 | 
                                       data_resettlement_oecd$year==2000 | data_resettlement_oecd$year==2019)
data_resettlement_3periods_t <- data_resettlement_3periods %>% dplyr::select(pol_openness_f, year)
table(data_resettlement_3periods_t$pol_openness_f, data_resettlement_3periods_t$year)

figure2a <- data_resettlement_3periods %>% 
  drop_na(pol_openness_f) %>% ggplot() + 
  geom_bar(aes(x=pol_openness_f, fill=factor(year), colour=factor(year)), position="dodge") +
  theme_bw() + xlab("") + ylab("Number of countries") + 
  scale_fill_manual("Year", values = c("grey90", "grey40", "gray0")) + 
  scale_colour_manual("Year", values = c("grey90", "grey40", "gray0")) + theme_minimal() +
  ggtitle("A. Resettlement policy") 

# Figure 2b
resettled_refugees <- data_resettlement_oecd %>% dplyr::group_by(year) %>% 
  dplyr::summarize(resettled = sum(reset, na.rm = TRUE))
resettled_refugees$year <- as.numeric(as.character(resettled_refugees$year))
resettled_refugees$resettled <- as.numeric(as.character(resettled_refugees$resettled))
figure2b  <- ggplot() + 
  geom_line(data=resettled_refugees, aes(y=resettled, x=year), color = "black", lwd=2) +
  theme_bw() + xlab("Year") + 
  scale_x_continuous(breaks=c(1980, 1990, 2000, 2010), limits=c(1980, 2019)) +
  scale_y_continuous(limits=c(0, 177000)) +
  ylab("Number of resettled refugees") + theme_minimal() +
  ggtitle("B. Resettlement admissions")

figure2 <- ggarrange(figure2a, figure2b, 
                     labels = c("", ""),
                     ncol = 2, nrow = 1, 
                     widths = c(1, 1.2))
jpeg("figure2.png", 
     width = 10, height = 5, units = "in", res=800)
plot(figure2)
dev.off()

# Figure 3
data_resettlement_oecd$AvgS_ExtCont_1 <- (data_resettlement_oecd$AvgS_ExtCont*(-1))+1
impic_compol <- data_resettlement_oecd %>% dplyr::filter(year<=2010) %>% dplyr::group_by(year) %>% 
dplyr::summarize(AvgS_ExtCont = mean(AvgS_ExtCont_1, na.rm = TRUE), 
                 Reset = mean(pol_openness, na.rm = TRUE))
impic_compol$year <- as.numeric(as.character(impic_compol$year))
impic_compol$AvgS_ExtCont <- as.numeric(as.character(impic_compol$AvgS_ExtCont))
impic_compol$Reset <- as.numeric(as.character(impic_compol$Reset))
blacktrans1 <- adjustcolor("black", alpha.f=0.3)
jpeg("figure3a.jpeg", 
     width = 8, height = 5, units = "in", res=800)
my_text <- "Control"
text1 = grid.text(my_text, x=0.94, y=0.37, gp=gpar(col="lightgray", fontsize=10, fontface="bold"))
my_text <- "Resettlement"
text2 = grid.text(my_text, x=0.94, y=0.67, gp=gpar(col="black", fontsize=10, fontface="bold"))
ggplot() + geom_jitter(aes(impic_compol$year, impic_compol$Reset), color="black", width=0.5, height=0.03) + geom_jitter(aes(impic_compol$year, impic_compol$AvgS_ExtCont), color="lightgray", width=0.5, height=0.03) + geom_hline(yintercept=c(0.25,0.5,0.75), linetype="dashed", color="gray") + ggtitle("A. Co-evolution of migration control and resettlement policy") + 
  geom_smooth(aes(impic_compol$year, impic_compol$Reset),method="loess", span=0.2, se=T, level=0.9, fill="black", color="black", na.rm=T, lwd=1.5)   + scale_x_continuous(name="Year", limits=c(1980,2013), breaks=c(1980,1990,2000,2010)) + scale_y_continuous(name="Policy openness", limits=c(0,0.75), breaks=c(0, 0.25,0.5,0.75,1)) +
  geom_smooth(aes(impic_compol$year, impic_compol$AvgS_ExtCont),method="loess", span=0.2, se=T, level=0.9,fill=blacktrans1, color="lightgray", na.rm=T, lwd=1.5) + theme_bw() + theme(plot.title=element_text(size=16, face="bold")) + annotation_custom(text1) + annotation_custom(text2) + theme_minimal()
dev.off()

resettlement_share <- data_resettlement_oecd %>% dplyr::filter(year<=2020) %>% dplyr::group_by(year) %>% 
  dplyr::summarize(resettlement = sum(reset, na.rm = TRUE), 
                   asylum = sum(asylum_seekers, na.rm = TRUE))
resettlement_share$share <- resettlement_share$resettlement/(resettlement_share$asylum+resettlement_share$resettlement)

data_resettlement_oecd$share <- data_resettlement_oecd$reset/(data_resettlement_oecd$asylum_seekers + data_resettlement_oecd$reset)
mix_ymean <- aggregate(data_resettlement_oecd$share, list(data_resettlement_oecd$year), mean, na.rm=T)
names(mix_ymean)[1] <- "year"
names(mix_ymean)[2] <- "share"
mix_ymean <- subset(mix_ymean, mix_ymean$share>=0 & mix_ymean$year<2020)
mix_ymean$ID1 <- NA
mix_ymean$ID1[mix_ymean$year==1980] <- 1
mix_ymean$ID1[mix_ymean$year==1985] <- 6
mix_ymean$ID1[mix_ymean$year==1990] <- 11
mix_ymean$ID1[mix_ymean$year==1995] <- 16
mix_ymean$ID1[mix_ymean$year==2000] <- 21
mix_ymean$ID1[mix_ymean$year==2005] <- 26
mix_ymean$ID1[mix_ymean$year==2010] <- 31
mix_ymean$ID1[mix_ymean$year==2015] <- 36
mix_ymean$ID1[mix_ymean$year==2019] <- 40

jpeg("figure3b.jpeg", 
     width = 8, height = 5, units = "in", res=800)
par (mar=c(3.5,4.5,2,1), mgp=c(2,.7,0), tck=-.01, las=1)   # for more TUFTE-STyle add also bty="l"
plot(resettlement_share$share, ylim=c(0,0.5), lwd=3, xaxt='n', yaxt='n', xlab="", ylab="", type="l",frame=F, cex=3)
axis(side = 1, at=mix_ymean$ID1, labels = mix_ymean$year, las=1, cex.axis=0.8) # at=41
axis(side = 2, at=c(0,0.1,0.2,0.3,0.4,0.5), labels = c("0%", "10%", "20%", "30%", "40%", "50%"), cex.axis=0.8, las=1) #
lines(mix_ymean$share, lwd=3, col="gray", add=T)
mtext("B. Resettlement share of humanitarian migration", side=3, adj=0, line=0.2, cex=1.1, font=1); 
mtext("Share of resettlement", side=2, adj=0.5, line=3.5, cex=0.9, las=3)
mtext("Year", side=1, adj=0.5, line=2.5, cex=0.9, font=1)
abline(h=c(0,0.1,0.2,0.3,0.4,0.5), lty=2, col="lightgray")
legend( x=30, y=0.39, bty = "n" ,
        legend=c("Total share","Average share"),
        col=c("black","gray"), lwd=3, pch=c(NA,NA) )
dev.off()

### Figure 4
data <- subset(data_resettlement_oecd, data_resettlement_oecd$share>=0 & data_resettlement_oecd$pol_openness>=0)
wth1 <- plm(reset_by_pop ~ as.factor(pol_openness), data=data, index=c("country", "year"), model="within")
btw1 <- plm(reset_by_pop ~ as.factor(pol_openness), data=data, index=c("country", "year"), model="between")
wth1a <- plm(share ~ as.factor(pol_openness), data=data, index=c("country", "year"), model="within")
btw1a <- plm(share ~ as.factor(pol_openness), data=data, index=c("country", "year"), model="between")
plot1 <- plot_models(wth1,btw1, show.values=T, title="A. Resettlement admissions (per 100,000 inh.)", colors=c("black", "grey50"),vline.color="grey90", line.size=1.2, show.p=T, axis.labels=c("ad hoc & quota", "quota","ad hoc"), legend.title = "Model", vcov.type = "HC3", m.labels = c("Within-country", "Between-country"))  + theme_blank()
plot2 <- plot_models(wth1a, btw1a, show.values=T, title="B. Resettlement-share of humanitarian arrivals", colors=c("black", "grey50"), vline.color="grey90", dot.size=3, line.size=1.2, show.p=T, axis.labels=c("ad hoc & quota", "quota","ad hoc"), legend.title = "Model", vcov.type = "HC3", m.labels = c("Within-country", "Between-country")) + theme_blank() + ylim(-0.3,0.54)
jpeg("figure4.jpeg", 
     width = 12, height = 5, units = "in", res=800)
ggarrange(plot1,plot2)
dev.off()

## Figure A1: Resettlement share of refugee population over time
share <- data_resettlement_oecd %>% dplyr::group_by(year) %>% 
  dplyr::summarize(refugees = mean(Refugees, na.rm = TRUE),
                   resettlement = sum(reset, na.rm = TRUE))
share$share <- share$resettlement/share$refugees   
jpeg("figure_a1.jpeg", 
     width = 8, height = 5, units = "in", res=800)
par (mar=c(3.8,5,0.4,0.4), mgp=c(2,.7,0), tck=-.01, las=1)
par(mfrow=c(1,1))
plot(share$year, share$share, type='l', ylim=c(0, 0.015),yaxt='n', xlab="", ylab="", frame=F, lwd=2)  
axis(side = 2, at=c(0,0.005,0.010,0.015), labels = c("0%", "0.5%", "1%", "1.5%"), las=1) #
mtext("Share of resettled refugees", side=2, adj=0.5, line=3.5, cex=1.3, las=3)
mtext("Year", side=1, adj=0.5, line=2.5, cex=1.3, font=1)
abline(h=c(0,0.005,0.01), lty=2, col="lightgray")
dev.off()

### Models

## base models
m1 <- lm(pol_openness ~ pop + gdp_lag + unemp_lag + lrs + rrpp_vote + refugees_lag + asylum_requests_lag + country, data = data_resettlement_oecd) 
coef_m1 <- coeftest(m1, cluster=data_resettlement_oecd$country)
m2 <- lm(pol_openness ~ pop + gdp_lag + unemp_lag + lrs + rrpp_vote + refugees_lag + asylum_requests_lag  + country, data = sample_europe)  
coef_m2 <- coeftest(m2, cluster=sample_europe$country)
m3 <- lm(reset_by_pop ~ pop + gdp_lag + unemp_lag + lrs + rrpp_vote + refugees_lag + asylum_requests_lag  + country, data = data_resettlement_oecd)  
coef_m3 <- coeftest(m3, cluster=data_resettlement_oecd$data_resettlement_oecd)
m4 <- lm(reset_by_pop ~ pop + gdp_lag + unemp_lag + lrs + rrpp_vote +  refugees_lag + asylum_requests_lag + country, data = sample_europe)  
coef_m4 <- coeftest(m4, cluster=sample_europe$country)

# Table 1
stargazer(m1, m2, m3, m4, title="Model estimates resettlement policy openness and resettlement admissions", 
          ci.level=0.95, star.cutoffs = c(0.05, 0.01, 0.001), 
          se = list(coef_m1[,2],coef_m2[,2],coef_m3[,2],coef_m4[,2]),
          column.sep.width = "10pt",
          no.space=TRUE,
          omit="country",
          dep.var.labels=c("Policy openness", "Resettlement admissions"), 
          #single.row = TRUE,
          omit.stat = c("f", "ser", "aic", "bic","adj.rsq"),
          covariate.labels = c("Population size", "GDP per capita", "Unemployment", "Left-wing cabinet", "Right-wing cabinet", "Radical-right vote share", 
                               "Refugees worldwide", "Asylum requests", "Constant"))

# logistic models with binary DV
logit1 <- glm(quota ~ pop + gdp_lag + unemp_lag + lrs + rrpp_vote + refugees_lag + asylum_requests_lag + country, data = data_resettlement_oecd,  family = binomial(link = "logit"))
coef_logit1 <- coeftest(logit1, cluster=data_resettlement_oecd$country)
logit2 <- glm(quota ~ pop + gdp_lag + unemp_lag + lrs + rrpp_vote + refugees_lag + asylum_requests_lag + country, data = sample_europe,  family = binomial(link = "logit"))
coef_logit2 <- coeftest(logit2, cluster=data_resettlement_oecd$country)

# Table A3
stargazer(logit1, logit2, title="Logistic models of resettlement policies", 
          ci.level=0.95, star.cutoffs = c(0.05, 0.01, 0.001), 
          se = list(coef_logit1[,2],coef_logit2[,2]),
          column.sep.width = "10pt",
          no.space=TRUE,
          omit="country",
          dep.var.labels=c("Resettlement policy (1=yes)"), 
          covariate.labels = c("Population", "GDP per capita", "Unemployment", "Left-wing cabinet", "Right-wing cabinet", "Radical-right vote share", 
                               "Refugees worldwide", "Asylum requests", "Constant"),
          omit.stat = c("f", "ser", "aic", "bic","adj.rsq"))

## Between-country models
m1b <- plm(pol_openness ~ pop + gdp_lag + unemp_lag + lrs + rrpp_vote + asylum_requests_lag, data = data_resettlement_oecd, index=c("country", "year"), model="between")  
coef_m1b <- coeftest(m1b, cluster=data_resettlement_oecd$country) # 
m2b <- plm(pol_openness ~ pop + gdp_lag + unemp_lag + lrs + rrpp_vote + asylum_requests_lag, data = sample_europe, index=c("country", "year"), model="between")  
coef_m2b <- coeftest(m2b, cluster=sample_europe$country) # 
m3b <- plm(reset_by_pop ~ pop + gdp_lag + unemp_lag + lrs + rrpp_vote + asylum_requests_lag, data = data_resettlement_oecd, index=c("country", "year"), model="between")  
coef_m3b <- coeftest(m3b, cluster=sample_europe$country) #
m4b <- plm(reset_by_pop ~ pop + gdp_lag + unemp_lag + lrs + rrpp_vote + asylum_requests_lag, data = sample_europe, index=c("country", "year"), model="between")  
coef_m4b <- coeftest(m4b, cluster=sample_europe$country) # 

# Table A4
stargazer(m1b, m2b, m3b, m4b, title="Between-country models", 
          ci.level=0.95, star.cutoffs = c(0.05, 0.01, 0.001), 
          se = list(coef_m1b[,2],coef_m2b[,2],coef_m3b[,2],coef_m4b[,2]),
          column.sep.width = "10pt",
          no.space=TRUE,
          omit="country",
          dep.var.labels=c("Policy openness", "Resettlement admissions"), 
          #single.row = TRUE,
          omit.stat = c("f", "ser", "aic", "bic","adj.rsq"),
          covariate.labels = c("Population size", "GDP per capita", "Unemployment", "Left-wing cabinet", "Right-wing cabinet", "Radical-right vote share", 
                               "Asylum requests", "Constant"))

### Analysing asylum policy mix
comb1 <- lm(pol_comb ~ rrpp_vote + asylum_requests_lag + country, data = data_resettlement_oecd)  
coef_comb1 <- coeftest(comb1, cluster=data_resettlement_oecd$country)
comb2 <- lm(pol_comb ~ rrpp_vote + asylum_requests_lag + country, data = sample_europe)  
coef_comb2 <- coeftest(comb2, cluster=sample_europe$country)

# Table A5
stargazer(comb1, comb2, title="Model estimates combined indicator (resettlement and border control policies)", 
          ci.level=0.95, star.cutoffs = c(0.05, 0.01, 0.001), 
          se = list(coef_comb1[,2],coef_comb1[,2]),
          column.sep.width = "10pt",
          no.space=TRUE,
          omit="country",
          dep.var.labels=c("Combined indicator"), 
          #single.row = TRUE,
          covariate.labels = c("Radical-right vote share", "Asylum requests", "Constant"),
          omit.stat = c("f", "ser", "aic", "bic","adj.rsq"))

# Figure 5
pol_comb <- subset(data_resettlement_oecd, data_resettlement_oecd$pol_comb>=-1)
lrs1 <- lm_robust(pol_comb ~ rrpp_vote + country, data = pol_comb)
coef_lrs1 <- coeftest(lrs1, cluster=data_resettlement_oecd$country)
lrs2 <- lm_robust(pol_comb ~ asylum_requests_lag + country, data = pol_comb)
coef_lrs2 <- coeftest(lrs2, cluster=data_resettlement_oecd$country)
jpeg("figure5.jpeg", 
     width = 11, height = 5.5, units = "in", res=800)
par (mar=c(3.8,5,0.4,0.4), mgp=c(2,.7,0), tck=-.01, las=1)
par(mfrow=c(1,2))
cplot(lrs1, "rrpp_vote", xlab = "Radical-right vote share", ylab ="Asylum policy mix", rug=F)
cplot(lrs2, "asylum_requests_lag", xlab = "Asylum requests", ylab ="Asylum policy mix", rug=F)
dev.off()

# Figure A2
data_resettlement_oecd$lrs <- factor(data_resettlement_oecd$lrs, levels = c("left", "center", "right"))
lrs1 <- lm_robust(pol_openness ~ lrs + country, data = data_resettlement_oecd)
lrs2 <- lm_robust(reset_by_pop ~ lrs + country, data = data_resettlement_oecd)  
jpeg("figureA2.jpeg", 
     width = 11, height = 5.5, units = "in", res=800)
p1 <- plot_model(lrs1, type = "eff", terms = "lrs") + ylim(0, 1) + xlab("Government ideology") + ylab("Policy openness") + theme_bw() + ggtitle("")
p2 <- plot_model(lrs2, type = "eff", terms = "lrs") + ylim(0, 13) + xlab("Government ideology") + ylab("Resettlement admissions") + theme_bw() + ggtitle("")
grid.arrange(p1, p2, nrow = 1)
dev.off()

# Figure A3
y2019 <- subset(data_resettlement_oecd, data_resettlement_oecd$year==2019)
gdp1 <- lm_robust(pol_openness ~ gdp_lag, data = y2019)  
gdp2 <- lm_robust(reset_by_pop ~ gdp_lag, data = y2019)  
jpeg("figureA3.png", 
     width = 11, height = 5.5, units = "in", res=800)
p1 <-plot_model(gdp1, type = "pred", terms="gdp_lag", show.values=T, show.data=T) + theme_bw() + xlab("GDP per capita (log-scale)") + ylab("Policy openness") + ggtitle("")
p2 <-plot_model(gdp2, type = "pred", terms="gdp_lag", show.values=T, show.data=T) + theme_bw() + xlab("GDP per capita (log-scale)") + ylab("Resettlement admissions") + ggtitle("")
grid.arrange(p1, p2, nrow = 1)
dev.off()

## correlations between resettlement openness and control restrictiveness
cor(data_resettlement_oecd$pol_openness, data_resettlement_oecd$AvgS_ExtCont, use="complete.obs")

func <- function(xx)
{
  return(data.frame(COR = cor(xx$pol_openness, xx$AvgS_ExtCont_1, use="complete.obs")))
}

ddply(data_resettlement_oecd, .(data_resettlement_oecd$country), func)


